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It is shown that homogeneous Rayleigh-Benard flow, i.e., Rayleigh-Benard turbulence with peri- 
odic boundary conditions in all directions and a volume forcing of the temperature field by a mean 
gradient, has a family of exact, exponentially growing, separable solutions of the full non-linear sys- 
tem of equations. These solutions are clearly manifest in numerical simulations above a computable 
critical value of the Rayleigh number. In our numerical simulations they are subject to secondary 
numerical noise and resolution dependent instabilities that limit their growth to produce statistically 
steady turbulent transport. 
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Much effort has been expended in recent decades in 
addressing the problem of heat transfer in Rayleigh- 
Benard thermal convection cells W^. An asymptotic 
high Rayleigh number heat transport scaling behavior 
Nu ~ Ra^/^ (perhaps with logarithmic modifications) 
has been conjectured as the ultimate regime 2, 3, 4, 5, 3- 
While current experimental data for high Rayleigh num- 
bers are controversial [sL Igl lTolllllfl^ . numerical simula- 
tions have not been very effective in studying this regime 
because of difficulties in dealing with the huge number 
of degrees of freedom engendered when Rayleigh num- 
bers reach the order of at least 10^^. Recently, some of 
us 1^, 1^ have studied a tri-periodic convective cell, or 
homogeneous Rayleigh-Benard (HRB) system, in order 
to bridge such difficulties, and to investigate the prop- 
erties of the convective cell once the effect of boundary 
layers has been eliminated. Model systems of this sort 
with hyperviscosity were first investigated computation- 
ally by Borue and Orszag , and later by Celani et al 
(in two spatial dimensions) ^14] with hyperviscosity and 
extra large scale dissipation as well. 

In this letter we point out some peculiar properties 
of the HRB model that are particularly striking in the 
low Rayleigh number regime, the opposite regime from 
that studied in refs 0, Q ■ First we display a family of 
exact, exponentially growing, separable solutions of the 
full non-linear HRB-system. We show that these solu- 
tions arc clearly manifest in direct numerical simulations 
in the Rayleigh number regime above a computable crit- 
ical value. Then by way of a careful numerical precision 
study we show that these may be robust and attracting 
solutions of the full system of partial differential equa- 
tions. 



Here we would like to anticipate that recently a clev- 
erly conducted series of experiments was designed such 
as to reduce the influence of top and bottom plates on 
the physical core of thermal convection 15j. In these 
experiments the temperature gradient in the bulk of the 
cell is not imposed but rather, as in fixed-flux convection 
[l6j |. measured as a dependent parameter. Interestingly, 
the Nu '-^ Ra}/'^ and Re ~ Ra}^^ scalings observed are 
consistent with HRB simulations in [s, 6J. 

The system to be studied is described in terms of the 
following partial differential equations 

ttt + It • Vu + Vp = j/Au + kag^, (1) 

AT 

Ot + u-ve = KAe + —u^, (2) 

H 

where u = {u^, Uy, u^) is an incompressible velocity 
field, V • It — 0, v and k and ag are respectively the 
kinematic viscosity, thermal diffusivity and the thermal 
expansion coefficient times the acceleration due to grav- 
ity. These equations are used to describe the evolution 
of the velocity field in a triply-periodic cubic volume 
[0, i/]'^ in the presence of a temperature field T{x,t) = 
T{x)-{-6{x, t). The temperature is expressed as a fluctua- 
tion 9 with respect to a mean profile T{x) that is imposed 
to be equal to the mean conductive temperature profile 
in such a Rayleigh-Benard cell; i.e. linear and of the form 
T{x) = -zAT/H + ^AT. 

When nondimensionalizing lengths with H, velocities 
with k/H, and temperatures with AT, equations Q and 
(PI can be rewritten as 

Ut + U -WU + VP ^ Pr(Ai7 + kRae), (3) 

Gi + ij-ve = Ae + w, (4) 
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where Q, P, and U = (u, w, w) are the dimensionless 
temperature, pressure, and velocity fields, respectively, 
and Pr = and Ra = agH^AT (vk) are the Prandtl 
and Rayleigh numbers. In this system periodic boundary 
conditions are imposed on all the dependent variables on 
the cube [0, 1]'^. We consider the equations of motion 
and 0) with spatially mean-zero initial data for all 
U and Q so that solutions subsequently remain spatially 
mean-zero at all times. 

With these boundary conditions there is a family of 
nonlinear "separable" solutions of (PJ and where U, 
Q and P are functions only of x, y and t but not of the 
vertical coordinate z. To see this, let v = m{x,y,t) + 
iv{x,y,t) and P = q{x,y,t). Then the divergence-free 
velocity v — (m, v), and q, w and O satisfy 

Vt+v - Vv + \7q = PrAv , (5) 
wt+v-\7w = Pr(Au; -I- RaO) , (6) 
Qt + v-VQ = AQ + w. (7) 

Equation © is the unforced two-dimensional Navier- 
Stokes equation whose solutions decay to zero exponen- 
tially in time 0|. As v decays away, equations ® and 
^ admit exact solutions of the form 

(S'(:;:;?))^(e:)^"»'"<'-+'»''+«' <«' 

with growth rate 

A = -i(Pr + l)k^ + iy(P7TT)2FT4R7(Ra~~fc4) (9) 

and where </> is an arbitrary phase. Therefore for Ra 
above a critical value Raf. — k'^ — (27r)^ the expo- 
nential solutions are unbounded; i.e. A > 0. Since 
k = 2tt {uxTriy), where n^x,y} are (±)-integer wave- 
numbers, we can re-label the exponent A as \(nx,ny). 
Degenerate values of A(nx,«j/) are possible, correspond- 
ing to different combinations of Ux, Uy. The number of 
positive A grows asymptotically Ra^^^. These solu- 
tions actually transport unlimited heat because Nu ^ 
(wQ) ^ exp(2At). Such runaway solutions and their 
possible instabilities may be actually the dominant fea- 
tures of the simulations in |0 , and in |^ where 
a Nu ~ Ra^^"^ dependence was observed. Therefore the 
relevance of computations on periodic domains as models 
for the bulk of systems with essential boundaries may be 
arguable. Indeed, Borue and Orszag [l^l seemed to hint 
at this issue when they remarked "It turns out that in ho- 
mogeneous convection [the heat transport, temperature 
and velocities] are strongly fluctuating and intermittent 
in time. This fact makes reliable measurements of [these 
variables] difficult." 

One of the purposes of this letter is to show that so- 
lutions in © appear to be attracting and may dominate 
the dynamics in numerical integrations of equations 
and Q on fully periodic domains at low Ra numbers. 
Resolution and precision difficulties, however, prevent us 
from drawing firm conclusions about the physics behind 



secondary instabilities that limit the growth of the run- 
away solutions to produce statistically steady turbulent 
transport in simulations at higher Ra. This remains an 
open question subject to further, deeper, investigations. 
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FIG. 1: Derivative with respect to time of the logarithm 
of Nu{t), in a homogeneous convective system at Ra = 
4.5 X 10"* ~ 30 • Rac. The DNS is implemented through a 
LBE algorithm on a cubic grid with resolution 240^, as in 
[^1^. Horizontal lines, from top to bottom, correspond re- 
spectively to values 2A(0, 1), 2A(1, 1), 2A(0,2), 2A(0,2) and 
the zero level. The time, here and in the following figures, is 
dimensionless as for the set of Eq. I^-Q, i.e., it has been nor- 
malized by the thermal difi^usion time across the box, /k. 

Indeed while unlimited heat transport could be ex- 
pected because of the solutions ©, even at moderately 
low Rayleigh numbers simulations of the homogeneous 
convective system display a statistically stationary be- 
havior where the growing modes, when they appear, 
break up due to rapid destabilization (we limit our in- 
vestigation to the Pr — 1 case). Fig. ^ shows the time 
derivative of the logarithm of the Nussclt number (i.e. 
the volume average Nu{t) = {wQ)v) for Ra ~ 4.5 x 10*. 
The growth rate of Nu{t) appears to bounce between 
some of the admissible exponential modes, although in 
these simulations the fastest growing modes, A(0, ±1) 
and A (±1,0), are never reached. The numerical results 
in Fig. n were obtained by means of a Lattice Boltzmann 
Equation (LBE) algorithm at a relatively high resolution, 
240'^ (see |g] for details). The value of Ra adopted here is 
a little below, but of the same order as, the lower values 
of Ra for data analyzed in . 

To better understand the relevant features of the dy- 
namics of the exponentially growing solutions we have 
performed a new series of direct numerical simulations 
(DNS) at values of Ra only slightly above the critical 
value {Ra > Rac) where just one distinct positive value 
of A exists. These integrations were performed by means 
of a fully de-aliased pseudo-spectral algorithm that al- 
lows for more flexibility. It allows the adjustment of the 
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time step size that is implicitly fixed by the spatial grid in 
the LBE, and gives clearer control of the scales involved 
in the dynamics. Furthermore, because we are interested 
here in the low-Ra regime, it is reasonable to perform 
numerical simulations with lower resolutions (i.e. 32'^ or 
64'^). Nevertheless we caution that in case of unlimited 
exponential growth any spatial resolution may be insuffi- 
cient at some point in time. Fig. |21 displays the temporal 









L/ 

\ 


v 




{/ 




i 






\ 
















i 









0.5 1 1.5 

time 



FIG. 2: Linear-Log plot of the spatial rms- value of the three 
velocity component and the thermal fluctuation in the DNS 
at resolution 32^, Wrms (solid), Urms (dashed), Vrms (dotted). 

behavior of the global (spatial) rms-values of the three 
velocity components. For clarity the temperature O has 
been omitted from the figure because it is strongly cor- 
related with w, almost coincident with it. As expected 
from the unstable analytic solution, w grows at an expo- 
nential rate oc exp (Xt) while the horizontal components 
(m, v) rapidly decrease. The plateaux in Fig. and the 
linear growth (on the lin-log scale) in Fig. |21 indicate the 
presence of exponentially growing solutions, but it is dif- 
ficult to determine the nature of the sudden departures 
from this state. Secondary instabilities are necessary to 
limit the growth of the runaway solutions and produce z- 
dcpcndcnt states that represent steady turbulent trans- 
port. However, Fig. |31 is a comparison between float- 
ing point and double precision calculations for the mean 
squared velocity components indicating the sensitivity of 
the exponential solutions to random perturbations gen- 
erated by round-off noise and discretization errors. 

The simulations indicated that the breakdown of the 
exponentially growing solutions for w and O and the ex- 
ponentially decaying solutions for u and v is first signaled 
in these horizontal components. It is not clear at this 
point what arc the relevant scales involved in this pro- 
cess for Ra > Ruc, although we observe that the growth 
of the horizontal components destabilizes the exponen- 
tially growing modes, leading to a fast redistribution of 
energy at all scales. At the peak, when w and Q reach 
their maxima, almost flat energy spectra are produced 
irrespective of the resolution adopted. Subsequently the 
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FIG. 3: A comparison between floating point (solid) and dou- 
ble precision (dashed) calculations for the root mean squared 
velocity component w, and {u^ +v'^y^^ versus time, with spa- 
tial resolution 32'^ and 2nd order Adams-Bashforth as time 
marching algorithm. 

high wave- vector dissipation takes over and the process is 
repeated, repeating the previous exponential growth. In 
summary, our simulations suggest that the exponential 
solutions may in fact be attracting for a broad class of 
initial conditions, although subject to some finite ampli- 
tude instabilities. 

There are other problems where extreme limits of mod- 
els introduce runaway solutions such as those appear- 
ing here. A distinct example is zero-Prandtl number 
Rayleigh-Benard convection p^ . IT^ where the singular 
limit of the Boussinesq equations admits exponentially 
growing solutions even in the presence of rigid (albeit 
free-slip) boundaries. It was long ago observed in double 
diffusive convection in the absence of rigid boundaries 
that linearly unstable modes representing 'salt fingers' 
are exact solutions of the nonlinear equations |20l |21| . 

Another model implemented with fully periodic condi- 
tions to avoid boundary layers, in which a similar nonlin- 
ear separation of variables appears, is shear-driven turbu- 
lence. In that model the fluctuations about an imposed 
mean shear flow iSy obey 

ut + u-Wu + Syda:U + iSv + Vp^ i^Au. (10) 

Fully periodic conditions cannot be implemented directly 
here, though, due to the presence of the incompatible op- 
erator Sydx (explicitely non-periodic in y). It was noted 
that periodic conditions can be imposed on indepen- 
dent variables x' = x — Syt, y' — y, z' — z and t' = t, and 
this transformation has been used to perform numerical 
simulations of "homogeneous shear flow" 23]. 

These periodic conditions also allow for an exact non- 
linear separation of the cross-stream and stream-wise 
components. Indeed, the change of variables implies 

W ^ \7' — jStdx' and Ut + SydxU — Ut' . (11) 
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Thus, acting on functions only of y' , z' and t' , the op- 
erators V = V = S^dy, + ka^/ and A = A' = 5^, + dl, . 
Hence the system pO|l separates for solutions depending 
only on y\ z' and t' : the two-dimensional divergence- free 
velocity fields v = jv{y' , z', t')+kw{y' , z' , t') and pressure 
p{y',z',t') satisfy the unforced Navier-Stokes equation 

Vt + v- Vv + Vp = vA'v, (12) 

and the stream-wise component u{y' , z' ,t') evolves ac- 
cording to the linear inhomogeneous equation 

ut + v ■V'u + Sv = I'A'u. (13) 

These equations do not (apparently) support unbounded 
exponentially growing fields but they do display non- 
normal transient growth among their fully nonlinear ex- 
act solutions. Indeed, the decaying solutions of (|12|l 

v = n e-'^('=2+'=3)t' sin(A:22;' + hz' + 0) [j/cs - k/ca] (14) 

produce a stream-wise flow of the form 

u={U- nSk2t') e-'^('=2+fe^)t' sm{k2y' + hz' + (f) (15) 

where U and are set by initial conditions. The peak 
amplitude Upeak — —e^^VlSkj,/ u{k2 -f fc|) (when U ~ Q) 
may be extremely large at high Reynolds number when 
Sk^/v{k2 + ^i) ^ 1- This is consistent with the be- 
havior reported by Pumir and Shraiman based on their 
direct numerical simulations |2,l| |: "The transient regime 
is characterized by a violent growth of the kinetic energy 
. . . While this growth eventually stops . . . the turbulent 
regime that follows exhibits large fluctuations of spatially 



averaged quantities . . . Because of the unusually large 
level of fluctuations, very long runs are necessary to get 
steady averages, which explains why we chose to work at 
moderate resolution." 

We conclude that the imposition of periodic bound- 
ary conditions may admit an exact nonlinear separability 
that allows for larger fluctuations than are possible in the 
presence of rigid boundaries. In full three-dimensional 
simulations of some of these systems, secondary instabil- 
ities are the only limiting processes that can lead to fi- 
nite statistically steady turbulent transport. Among the 
delicate — and currently unresolved — issues is the ques- 
tion of how sensitive the statistics of high Ra (or Re) 
simulations may be to numerical discretization and noise. 
Nevertheless the HRB model, although physically un- 
realizable because of the boundary conditions, has stim- 
ulated new and interesting experiments where the ef- 
feet of the thermal boundaries has been reduced to reveal 
the Nu ~ Ra}/^ scaling observed in the high Rayleigh 
number simulations 5, &\. 
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